NASA / TM-1999-209365 
ARL-TR-2012 



A Method for Calculating Strain Energy 
Release Rates in Preliminary Design of 
Composite Skin / Stringer Debonding 
Under Multi- Axial Loading 


Ronald Krueger 

National Research Council 

Langley Research Center, Hampton, Virginia 


Pierre J. Minguet 

The Boeing Company, Philadelphia, Pennsylvania 


T. Kevin O'Brien 

U. S. Army Research Laboratory 
Vehicle Technology Directorate 

Langley Research Center, Hampton, Virginia 


July 1999 



The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated to 
the advancement of aeronautics and space 
science. The NASA Scientific and Technical 
Information (STI) Program Office plays a key 
part in helping NASA maintain this important 
role. 

The NASA STI Program Office is operated by 
Langley Research Center, the lead center for 
NASA's scientific and technical information. The 
NASA STI Program Office provides access to the 
NASA STI Database, the largest collection of 
aeronautical and space science STI in the world. 
The Program Office is also NASA's institutional 
mechanism for disseminating the results of its 
research and development activities. These 
results are published by NASA in the NASA STI 
Report Series, which includes the following 
report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant 
phase of research that present the results of 
NASA programs and include extensive 
data or theoretical analysis. Includes 
compilations of significant scientific and 
technical data and information deemed to 
be of continuing reference value. NASA 
counterpart of peer-reviewed formal 
professional papers, but having less 
stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary 
or of specialized interest, e.g., quick release 
reports, working papers, and 
bibliographies that contain minimal 
annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or co-sponsored by 
NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific 
and technical material pertinent to NASA's 
mission. 

Specialized services that complement the STI 
Program Office's diverse offerings include 
creating custom thesauri, building customized 
databases, organizing and publishing research 
results ... even providing videos. 

For more information about the NASA STI 
Program Office, see the following: 

• Access the NASA STI Program Home Page 
at http://www.sti.nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI Help 
Desk at (301)621-0134 

• Phone the NASA STI Help Desk at 
(301) 621-0390 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7121 Standard Drive 

Hanover, MD 21076-1320 


NASA / TM-1999-209365 
ARL-TR-2012 



A Method for Calculating Strain Energy 
Release Rates in Preliminary Design of 
Composite Skin / Stringer Debonding 
Under Multi- Axial Loading 


Ronald Krueger 

National Research Council 

Langley Research Center , Hampton, Virginia 

Pierre J. Minguet 

The Boeing Company, Philadelphia, Pennsylvania 

T. Kevin O'Brien 

U. S. Army Research Laboratory 
Vehicle Technology Directorate 

Langley Research Center, Hampton, Virginia 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


July 1999 



Available from: 


NASA Center for AeroSpace Information (C ASI) 
7121 Standard Drive 
Hanover, MD 21076- 1 320 
(301)621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22 1 6 1 -2 1 7 1 
(703) 605-6000 



A METHOD FOR CALCULATING STRAIN ENERGY RELEASE 
RATES IN PRELIMINARY DESIGN OF COMPOSITE 
SKIN/STRINGER DEBONDING UNDER MULTI-AXIAL 

LOADING 


Ronald Krueger 1 , Pierre J. Minguet 3 , and T. Kevin O’Brien 2 
1 National Research Council Research Associate 
2 U.S. Array Research Laboratory, Vehicle Technology Directorate 
NASA Langley Research Center 
Hampton, VA 23681 
3 Boeing 

Philadelphia, PA 19142 


ABSTRACT 

Three simple procedures were developed to determine strain energy release rates, G, in 
composite skin/stringer specimens for various combinations of uniaxial and biaxial 
(in-plane/out-of-plane) loading conditions. These procedures may be used for parametric design 
studies in such a way that only a few finite element computations will be necessary for a study of 
many load combinations. The results were compared with mixed mode strain energy release rates 
calculated directly from nonlinear two-dimensional plane-strain finite element analyses using the 
vir tual crack closure technique. The first procedure involved solving three unknown parameters 
needed to determine the energy release rates. Good agreement was obtained when the external 
loads were used in the expression derived. This superposition technique, however, was only 
applicable if the structure exhibits a linear load/deflection behavior. Consequently, a second 
modified technique was derived which was applicable in the case of nonlinear load/deformation 
behavior. The technique, however, involved calculating six unknown parameters from a set of six 
simultaneous linear equations with data from six nonlinear analyses to determine the energy release 
rates. This procedure was not time efficient, and hence, less appealing. 

Finally, a third procedure was developed to calculate mixed mode energy release rates as a 
function of delamination lengths. This procedure required only one nonlinear finite element 
analysis of the specimen with a single delamination length to obtain a reference solution for the 
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energy release rates and the scale factors. The delamination was subsequently extended in three 
separate linear models of the local area in the vicinity of the delamination subjected to unit loads to 
obtain the distribution of G with delamination lengths. This set of sub-problems was solved using 
linear finite element analyses, which resulted in a considerable reduction in CPU time compared to 
a series of nonlinear analyses. Although additional modeling effort is required to create the local 
sub-model, this superposition technique is very efficient for large parametric studies, which may 
occur during preliminary design where multiple load combinations must be considered. 
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Composite materials, fracture mechanics, energy release rate, finite element analysis, virtual crack 
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INTRODUCTION 

Carbon epoxy composite structures are widely used by today's aircraft manufacturers to 
reduce weight. Many composite components in aerospace structures consist of flat or curved 
panels with co-cured or adhesively bonded frames and stiffeners. Testing of stiffened panels 
designed for pressurized aircraft fuselage has shown that bond failure at the tip of the frame flange 
is an important and very likely failure mode [1]. Comparatively simple simulation specimens 
consisting of a stringer bonded onto a skin were developed and it was shown in experiments that 
the failure initiated at the tip of the flange, identical to the failure observed in the full-size panels 
and frame pull-off specimens [2-7]. 

The overall objective of the current work is to develop a simple procedure to calculate the 
strain energy release rate for delaminations originating from matrix cracks in these skin/stringer 
simulation coupons for arbitrary load combinations. The total strain energy release rate would then 
be compared to critical values obtained from an existing mixed-mode failure criterion to predict 
delamination onset. This procedure could then be used for parametric design studies in such a way 
that only a few finite element computations would be necessary to evaluate bonded joint response 
due to many load combinations. Since energy is a quadratic function of the applied loads, simple 
superposition to add the energy release rates from separate load cases is not valid. Therefore, a 
simple quadratic expression is developed to calculate the strain energy release rate for any 
combination of loads [4], To validate this approach, results obtained from the quadratic expression 
are compared to mode I and mode II strain energy release rate components, which are calculated 
from nonlinear two-dimensional plane-strain finite element analyses using the virtual crack closure 
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technique [8, 9], 

Three simple procedures are developed to determine strain energy release rates, G, in 
composite skin/stringer specimens for various combinations of uniaxial and biaxial 
(in-plane/out-of-plane) loading conditions. The first procedure involved solving three unknown 
parameters needed to determine the energy release rates. This superposition technique, however, 
was only applicable if the structure exhibits a linear load/deflection behavior. Consequently, a 
second modified technique is derived which is applicable in the case of nonlinear load/deformation 
behavior. A third procedure is developed to calculate mixed mode energy release rate as a function 
of delamination length. This procedure requires only one nonlinear finite element analysis of the 
specimen with a single delamination length to obtain a reference solution for the energy release 
rates and the scale factors. 


BACKGROUND 

Previous investigations of the failure of secondary bonded structures focused on loading 
conditions as typically experienced by aircraft crown fuselage panels. Tests were conducted with 
specimens cut from a full-size panel to verify the integrity of the bondline between the skin and the 
flange or frame [1], However, these panels were rather expensive to produce and there is a need 
for a test configuration that would allow detailed observations of the failure mechanism at the 
skin/flange interface. A simpler specimen configuration was proposed in reference 2. The 
investigations focused on the failure mechanisms of a bonded skin/flange coupon configuration 
loaded in bending [2-5]. In many cases, however, composite structures may experience both 
bending and membrane loads during in-flight service. Damage mechanisms in composite bonded 
skin/stringer structures under monotonic tension, three-point bending, and combined 
tension/bending loading conditions were investigated in references 6 and 7. An analytical 
methodology was also developed to predict the location and orientation of the first transverse 
matrix crack based on the principal transverse tension stress distribution in the off axis plies nearest 
the bondline in the vicinity of flange tip. The prediction of delamination onset was based on energy 
release rate calculations. 

The specimens tested in references 6 and 7 consisted of a bonded skin and flange assembly 
as shown in Figure 1. Both the skin and the flange laminates had a multidirectional lay-up made 
from IM6/3501-6 graphite/epoxy prepreg tape with a nominal ply thickness of h =0. 188 mm. The 
skin lay-up, consisting of 14 plies, was [0/45/90/-45/45/-45/0] s and the flange lay-up, consisting 
of 10 plies, was [45/90/-45/0/90] s . The measured bondline thickness averaged 0.102 mm. 
Specimens were 25.4-mm wide and 203.2-mm long. Typical material properties for the composite 


3 



tape and the adhesive material used in the analysis were taken from reference 2 and are summarized 
in Table 1. 

The specimens were subjected to pure tension, three-point bending, and combined axial 
tension and bending loads. A schematic of the deformed specimen geometries, the boundary 
conditions, and the loads corresponding to the first damage observed are shown in Figure 2. In 
the combined axial tension and bending load case, a constant axial load, P, was applied in a first 
load step while transverse loads remained zero. In a second load step, the axial load was kept 
constant while the load orientation rotated with the specimen as it deformed under the transverse 
load. The tests were terminated when the flange debonded unstably from one of the flange tips. 
Damage was documented from photographs of the polished specimen edges at each of the four 
flange comers identified in Figure 3(a). Typical damage patterns, which were similar for all three 
loading configurations, are shown in Figure 3(b) and (c). Comers 1 and 4 and comers 2 and 3 
had identical damage patterns. At comers 1 and 4, a delamination running in the 90745° flange ply 
interface (delamination A) initiated from a matrix crack in the 90° flange ply as shown in 
Figure 3(b). At longer delamination lengths, new matrix cracks formed and branched into both the 
45° ply below the delaminated interface as well as the 90° flange ply above the interface. At comers 
2 and 3 a matrix crack formed at the flange tip in the 90° flange ply that subsequently ran through 
the lower 45° flange ply and the bondline into the skin as shown in Figure 3(c). Subsequently, a 
split (delamination Bl) formed from the tip of that matrix crack within the top 0° skin ply and in 
some cases, a second delamination (delamination B2) was observed below the first in the top 
0745° skin ply interface. 

In previous investigations, stress analyses were used to predict the location and orientation 
of the first transverse matrix crack based on the principal transverse tension stress distribution in 
the off axis plies nearest the bondline in the vicinity of the flange tip [6,7]. A comparison of the 
trajectories of the maximum principle tension stress with the damage patterns shown in Figures 
3(b) and (c) indicated that the matrix crack starts to grow perpendicular to the trajectories. For all 
three loading conditions, maximum principal tensile stresses in the 90° ply closest to the bondline, 
computed for applied loads at damage onset, were almost identical and exceeded the transverse 
tension strength of the material. Subsequent finite element analyses of delamination growth from 
these matrix cracks were performed using the virtual crack closure technique. However, because 
the specimen geometry and loadings required nonlinear analyses, this was a computationally 
intensive process. 
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ANALYSIS FORMULATION 


FINITE ELEMENT MODEL 

In the current investigation the finite element (FE) method was used to analyze the test 
specimens for each loading case. The goal of this analysis is to evaluate strain energy release rate 
components at the delamination tip using the virtual crack closure technique [8,9], To develop a 
simple procedure to calculate the strain energy release for delaminations originating from matrix 
cracks, it was reasonable to focus only on one damage pattern during die investigation. Therefore, 
only a FE model of a specimen with a delamination running in the 90745° flange ply interface, 
corresponding to Figure 3b, was developed and loads and boundary conditions were applied to 
simulate the three load cases. The two-dimensional cross section of the specimens was modeled 
using quadratic eight-noded quadrilateral plane strain elements (see Figure 4) and a reduced (2x2) 
integration scheme was used for these elements. For the entire investigation, the ABAQUS finite 
element software was used [10]. 

An outline and two detailed views of the FE model are shown in Figure 4. A refined mesh 
was used in the critical area of the 90° flange ply where matrix cracks and delaminations were 
observed in the test specimens. Outside the refined mesh area, all plies were modeled with one 
element through the ply thickness. Two elements were used per ply thickness in the refined region, 
except for the first three individual flange plies above the bondline and the skin ply below the 
bondline, which were modeled with four elements. Three elements through-the-thickness were 
used for the adhesive film. Based upon the experimental observations shown in Figure 3b, the 
model included a discrete matrix crack and a delamination. The initial matrix crack was modeled 
perpendicular to the flange taper, as suggested by the microscopic investigation as well as the 
stress analysis, which showed that the matrix crack starts to grow perpendicular to the trajectory of 
the maximum principle tension stress [6,7], Damage was modeled at one flange tip as shown in 
Figure 4. The mesh used to model the undamaged specimen, as discussed in reference 6 and 7 , 
was employed at the opposite taper. The model consisted of 6977 elements and 21486 nodes and 
had 42931 degrees of freedom. 

For the combined tension and bending load case, performed in NASA Langley's axial 
tension and bending test frame [11,12], the top grip, the load cell, and the load pin were modeled 
using three-noded quadratic beam elements as shown in Figures 2c and 5, to accurately simulate 
the combined tension and bending loads applied [6,7]. The beams were connected to the 
two-dimensional plane strain model of the specimen using multi-point constraints to enforce 
appropriate translations and rotations. As shown in Figure 5, nodes 1-29 along the edge of the 
plane strain model (x =101.6 mm) were constrained to move as a plane with the same rotation as 
beam node A. To be consistent with the actual tests, a constant axial load, P, was applied in a first 
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load step while transverse loads remained zero. In a second load step, the axial load was kept 
constant while the load orientation rotated with the specimen as it deformed under the transverse 
load. During the tests, the maximum specimen deflections under the transverse load were recorded 
at the top grip contact point. In the FE simulation a prescribed displacement, v, was applied which 
corresponded to the recorded transverse stroke. For the beam model of the steel parts (top grip, 
load cell, and load pin), a Young’s Modulus of 210 GPa and a Poisson's Ratio of 0.3 were used 
as material input data. A rectangular beam cross section was selected to model the square cross 
section of the top grip (/ =1.87 x 10 6 mm 4 ) and load pin (/ =1.4 x 10 6 mm 4 ) and a circular beam 
cross section was used to model the cylindrical load cell (/ =8.37 x 10 3 mm 4 ). 

When applying two dimensional plane strain FE models it is assumed that the geometry, 
boundary conditions and other properties are constant across the entire width of the specimen. The 
current model, thus, may not always capture the true nature of the problem. As shown in 
Figure 3, the delamination pattern changed from comer 3 to comer 4 from a delamination 
running in the 90745° interface to a delamination propagating between the adhesive film and the top 
0° ply of the skin. This is a three dimensional effect and can not be accounted for in the current 
plane strain model. 


VIRTUAL CRACK CLOSURE TECHNIQUE 

The Virtual Crack Closure Technique (VCCT) described in references 8 and 9 was used to 
calculate strain energy release rates for the delaminations. The mode I and mode II components of 
the strain energy release rate, G t and G n , were calculated as (see Figure 6) 

and 

G »=-2^hK-“V) +x, i(“'<-7-)] ® 

where to is the length of the elements at the delamination tip, X' x and F.' are the forces at the 
delamination tip at node i, and « m ’ and v m ’ are the relative displacements at the corresponding node 
m behind the delamination tip as shown in Figure 6. S imilar definitions are applicable for the 
forces at node j and displacements at node L For geometrically nonlinear analysis, both forces and 
displacements were transformed into a local coordinate system (*', y'), that defined the normal and 
tangential coordinate directions at the delamination tip in the deformed configuration. The mode III 
component is identically zero for the plane strain case. Therefore, the total strain energy release 
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rate, G T , was obtained by summing the individual mode components as 

G J =G l +G n . (3) 

The data required to perform the VCCT in equations (1) to (3) were accessed directly from 
the ABAQUS binary result file to get better accuracy. The calculations were performed in a 
separate post processing step using nodal displacements and nodal forces at the local elements in 
the vicinity of the delamination front 

Care must be exercised in interpreting the values for G, and G n obtained using the virtual 
crack closure technique for interfacial delaminations between two orthotropic solids [13,14], For 
the current investigation, the element length Aa was chosen to be about 1/4 of the ply thickness, h, 
for the delamination in the 90°/45° flange ply interface. Note that for the FE model shown in 
Figure 4 Aa/h =0.181 for the element behind and Aa/h =0.25 for the element in front of the 
delamination tip. Therefore, the technique suggested in reference 8 was used to estimate the forces 
Xj' and Y t ' for the case of unequal element lengths at the delamination tip. For the further 
delamination growth a value of Aa/h =0.25 was used. 


ANALYTICAL INVESTIGATION 


SUPERPOSITION TECHNIQUE FOR LINEAR DEFORMATION BEHAVIOR 


The schematics of the specimen, boundary conditions, and three load cases (tension, 
bending and combined tension and bending) considered in this part of the study are shown in 
Figure 7. These boundary conditions and loads, however, do not represent the conditions applied 
during the experiments as given in Figure 2 of the previous section. This new set of boundary 
conditions was chosen to simplify the derivation of the superposition technique for linear 
deformation behavior. It was postulated that the specimen exhibits a linear load deflection behavior 
for the three load cases shown. Only linear finite element analyses were used. The boundary 
conditions applied were the same for all load cases. 

For a specimen subjected to a pure tension load P as shown in Figure 7(a), the energy 
release rate G p at the delamination tip can be calculated as 


G P 


EL E£p. 

2 ' BA 


(4) 


where C p is the compliance of the specimen and dA is the increase in surface area corresponding to 
an incremental increase in load or displacement at fracture [15]. For a specimen subjected to a 
bending load Q, as shown in Figure 7(b), the energy release rate G Q at the delamination tip can be 
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calculated accordingly as 



Q 2 d c g 

2 dA ' 


(5) 


If the external load, R, applied in the linear analysis is simply a fraction or multiple of the tension 
load P,R = nP, or the bending load Q, R = mQ, the energy release rate G R for the new load case 
may be obtained from the known values using 


G r — n 2 G P or G R — m 2 Gg. (6) 

In the case of a combined tension/bending load case as shown in Figure 7(c), where the external 
load is a combination of a fraction or multiple n of the tension load P and a different fraction or 
multiple m of the bending load Q,R = nP + mQ, we obtain 


QiP + mg) 2 3Q _ (n 2 / >2 +2mnP6+m 2 Q 2 j dCp 
* " 2 ’ dA ~ 2 ' dA' 


(7) 


Note that for a tension load, P, only. 


dC R _ dC f 


and for a bending load, Q, only. 


dC D dC r 


dA dA 

For the combined load case equation (7) can then be approximated by 


dA 


JL — 

dA 


G r = 


n 2 P 2 dCp 
2 ' dA 


+ 2mn 


PQdCt 
2 dA 


+ 


nj*Q^ 

2 dA ’ 


( 8 ) 


Using equations (4) and (5) yields 


G r = n 2 G P + 2mn ■ 


PQ 

2 


dC, 


dA 


— + m 2 G, 


Q' 


r PQ 


(9) 


where Gpq is a coupling term which has the dimension of an energy release rate. 

First, linear FE analyses of a simple tension and simple bending case are performed using 
VCCT to determine G,, G„ and G T . This allows calculation of the G P and G Q parameters in 
equation (9) for total G, and the G, and G n components. Then a single linear FE analysis of a 
combined tension and bending load case is performed using VCCT to obtain the G R parameter in 
equation (9) for G„ G„ and G T . Once these parameters are determined, then G PQ may be calculated 
for G,, G„ and G r . The parameters Gp, G Q and Gpq may now be used to calculate G R for G /; G„ and 
G T for other tension and bending load combinations. 

Mode I and mode II values were computed using VCCT for a delamination running in the 
90745° flange ply interface with a length equal to the length of the first element (a/h = 0.181) as 
shown in Figure 4. For the pure tension and bending loads shown in Figures 7(a) and (b), energy 
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release rates were also calculated using the analytical expressions of equation (6). In the example 
shown in Figure 8 for the tension load case, the parameter G p in equation (6) was computed for P= 
5.5 kN. The total energy release rate G T computed using VCCT and the superposed results are 
identical, since equation (6) is an exact closed form solution. Minor differences for the individual 
modes, that cannot be explained, are observed. For all permutations of P and Q loads, as shown in 
Figure 7(c), energy release rates for the combined load case were calculated using equation (9). In 
this investigation the parameter G p in equation (9) was calculated for a tension load P= 5.5 kN, G c 
was determined for a bending load Q= 1 12.5 kN and G PQ was obtained from one analysis of the 
combined tension and bending load. Energy release rates obtained from equation (9) were 
compared to mode I and mode II values calculated using VCCT as shown in Figure 9 for the case 
where a tension load P= 11.0 kN was applied and Q was varied. For the other permutations of 
loads the comparisons of only the total energy release rates, G ^ are shown in Figure 10. The good 
agreement of results confirms that the superposition technique derived in equation (9) is applicable, 
in combination with linear finite element analysis and VCCT to determine the unknown parameters, 
provided the structure shows a linear load/deflection behavior. 


A MODIFIED TECHNIQUE FOR NONLINEAR DEFORMATION BEHAVIOR 

For the investigation of the combined axial tension and bending load case as shown in 
Figures 2(c) and 5, nonlinear finite analyses were used since this allowed the axial load to rotate 
with the specimen as it deformed under the transverse load and accounted for the membrane 
stiffening effect caused by the axial load. In this case the superposition technique derived for the 
linear case in the previous section (equations (8) and (9)) is no longer applicable and a modified 
method needs to be developed. 

An analytical expression was suggested in reference 4 that is primarily a modification of 
equation (8) derived in the previous section. The external tension load, P, and bending load, Q, in 
the analytical expression were replaced with the local force resultant AA^and moment resultant 
yielding 

G = G^Ml + 2G mn M xx N xx +G nn A£, (10) 

where G^ and G M are unknown parameters determined from a pure tension and a pure bending 
load case and G^ is an unknown combined tension and bending parameter. The local force and 
moment resultants are calculated at the flange tip as shown in Figure 11. Location and calculation 
of the force and moment resultants. For improved accuracy, the terms related to the transverse 
shear force resultant, Q xy , were also included in expression (10) yielding 


9 




Unlike the linear case where a pure tension or a pure bending load case alone may be used to 
determine one of the unknown parameters, nonlinear analysis of the pure tension and pure bending 
load case yielded a combination of M a and at the flange tip due to the load eccentricity (tension 
load) and large displacements (bending load). Therefore, the constants Gy (i,j=m,n,q) could not be 
determined simply from the pure tension and bending load cases. Consequently, all six constants 
were calculated from a set of six simultaneous linear equations corresponding to six unique loading 
combinations solved previously, using nonlinear FE analyses. This yields G t ( k =1,...,6). 


"^1 I'M, 2 2 M X N X N 2 2 M X Q X 2N X Q X Qf \ TG^' 

G 2 m\ 2M 2 N 2 Nl 2 M 2 Q 2 2 M z Q 2 Ql G m 

G 3 _ Ml 2 M 3 N 3 Nl 2M 3 Q 3 2M 3 Q 3 Q\ G^ 

G 4 “ Ml 2M 4 N 4 Nl 2M 4 Q 4 2M 4 Q 4 Q 2 4 G^ 

G 5 Ml 2 M 5 N 5 Nl 2M 5 Q 5 2 M s Q 5 Ql G^ 

G 6 J |_M 2 2 M 6 N 6 Nl 2M 6 Q 6 2M 6 Q 6 Q 2 ] [_G qq _ 


(13) 


Further, the local force and moment resultants N xx , M xx , and Q xy for all six unique loading 
combinations were calculated at the flange tip using die equations shown in Figure 11 by 
integrating stresses determined in the nonlinear FE analyses yielding N v A/ k , and Q^{k = 

The system of six equations was then solved for the unknown G X] values. With the constants G Xj 
known, G could then be calculated from the force and moment resultants N xx , and M xx for any 
combined tension/bending load case using the technique described by equation (11). The term G is 
used here for the total energy release rate or for a mixed mode energy release rate component. 
Hence, the calculation of each of the individual modes G p G„ or G r requires a unique set of G t) 
constants each. This means that equation (13) needs to be solved individually for each fracture 
mode (I,II) before equation (1 1) is used to obtain the individual modes G p G u or G r 

The analytical expressions (10) and (11) were derived with the objective of developing a 
simple procedure to calculate the strain energy release rate if the specimen shows a nonlinear 
load/deflection behavior. The expressions may also be used if the specimen exhibits a linear 
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load/deflection behavior. Calculating the force and moment resultants and solving equation (13) to 
obtain a unique set of constants G-- for each fracture mode, however, appears to be cumbersome in 
this case because FE analysis needs to be performed for six unique combined load cases to 
determine the unknown parameters G y . In contrast, the use of expression (8) is simpler, because 
the external loads are known and only three load cases need to be analyzed to determine G p , G Q 
and G ' pq • 

The matrix equation (13), which contains the terms of local force and moment resultants 
N t , M v and <2 k , may become singular. For linear load/deflection behavior this will occur if at least 
one of the six load cases selected to calculate N t , Af k , and G k is not independent from the other 
cases, but simply a linear combination of any of them. For nonlinear load/deflection behavior it is 
not easily predictable under which circumstances the matrix might become singular. In both cases, 
however, six unique load cases need to be selected to avoid matrix singularity and solve equation 
(13) for the unknown parameters. 

The energy release rates were calculated using die modified method (equation (11)) for all 
permutations of axial loads, P, and transverse displacements, v max , shown in Figure 5. The 
unknown parameters Gy in equation (13) were obtained from nonlinear finite element analyses of 
six different unique load cases (P t = 0.0, v,= 30.9 mm; P 2 = 4.5 kN, v 2 = 7.5 mm; P 3 - 4.5 kN, v 3 = 
30.9 mm; P 4 = 9.0 kN, v 4 = 7.5 mm; P 5 = 9.0 kN, v 5 = 30.9 mm; P 6 =17.8 kN, v 6 = 30.9 mm). 
Calculated mixed-mode results were compared with the energy release rates obtained directly from 
nonlinear finite element analyses using VCCT as shown in Figure 12 for a case where only one 
axial load of P = 4.5 kN and multiple transverse displacements, v maT , were applied. As expected, 
the results were identical for the two cases which had been selected to determine the unknown 
parameters G iy For the other load combinations, G,, G„ and G r were in excellent agreement Total 
energy release rates calculated for all axial load and transverse displacement permutations are 
shown in Figure 13. For the remaining load combinations, calculated strain energy release rates 
differed by less than 5% when compared to results computed directly from nonlinear finite element 
analysis using VCCT. Good results, however, were only obtained if the six unique load 
combinations to determine the unknown parameters G ;j include the upper and lower limits of load 
combinations as shown in Figure 13. The modified method should be used to interpolate results 
for different load combinations. Extrapolation may lead to inaccurate results. 

Hence, it was possible to derive a technique which was applicable for nonlinear 
deformation of the specimen. The expression derived for die linear case was modified such that 
terms of the external forces were replaced by internal force and moment resultants. The energy 
release rates calculated using this technique seemed sufficiently accurate for preliminary design 
studies. However, while external forces are known, force and moment resultants at the flange tip 
need to be calculated analytically or computed from finite element analysis. For the current study of 
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the combined axial tension and bending load case, nonlinear finite analyses were used to c alculate 
the force and moment resultants. This requires about the same computational effort as directly 
computing the energy release rates from nonlinear analyses using the virtual crack closure 
technique. An additional effort is required to obtain the unknown parameters G Vj . The use of the 
technique as given in equation (11) may therefore become time consuming and less appealing for 
quickly calculating energy release rates for a large number of new load combinations from a set of 
known results. Furthermore, this process may have to be repeated for each new delamination 
length modeled to obtain the distribution of G p G n and G T as a function of delamination length. 
Consequently, another approach was developed for the simulation of delamination growth. 


SIMULATION OF DELAMINATION GROWTH 

The techniques developed in the previous sections focused on simple procedures to 
calculate the strain energy release rate for various combinations of loads from results previously 
computed for other load cases. A related problem is the simulation of delamination growth where 
mixed mode energy release rates need to be calculated as a function of delamination length, a. The 
shape of the G versus a curves for G p G u and G T yield information about stability of delamination 
growth and often dictate how these energy release rates are used to predict the onset of 
delamination [16]. During the nonlinear finite element analyses, the delaminations are extended and 
strain energy release rates are computed at virtual delamination lengths using the virtual crack 
closure technique. For preliminary design studies with several load cases of interest, delamination 
positions and lengths need to be checked continuously. Hence, the amount of computation time 
necessary may become excessive. Therefore fast and accurate alternatives need to be developed. 


REVIEW OF SIMULATED DELAMINATION PROPAGATION USING A SERIES OF 
NONUNEAR FINITE ELEMENT ANALYSES 

The schematics of the deformed geometries, the boundary conditions, and the loads 
examined in this part of the study are shown in Figure 2 for all three load cases. The boundary 
conditions considered in the simulations were chosen to model the actual test from references 6 and 
7 as closely as possible. For the tension and bending case, the mean loads reported for the point of 
damage initiation were applied. At this point, matrix cracks are likely to form. To be consistent 
with the combined axial tension and bending tests, a constant axial load, P = 17.8 kN, was applied 
in a first load step while transverse loads remained zero. In a second load step, the axial load was 
kept constant while the load orientation rotated with the specimen as it deformed under the 
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transverse load. In the FE simulation, a prescribed displacement was applied which corresponded 
to the average of the transverse stroke ( v = 31 mm) for which flange debond occurred [6,7]. 

The initial matrix crack was modeled on one flange tip perpendicular to the flange taper as 
suggested by the microscopic investigation and shown in Figure 3. The model of the discrete 
matrix crack and delamination is shown in Figure 4. During the nonlinear finite element analyses, 
the delaminations were extended and strain energy release rate components were computed as a 
function of delamination length using the virtual crack closure technique. The delamination lengths, 
a, were measured from the end of the initial matrix crack as shown in Figure 4. The delamination 
was extended in twelve increments up to about 0.6 mm (a/h = 3.2) which corresponds to a length 
where matrix crack branches were observed in the experiments as shown in Figure 3(b). The 
simulated delamination propagation therefore required 12 nonlinear FE analyses for each load case, 
consequently 36 analyses for all three load cases. The results plotted in Figures 14 through 16 
show that G n increases monotonically for all load cases while G t begins to level off at the longest 
delamination lengths [6,7]. These results were intended as reference solutions to be compared with 
results from the superposition method in the following section. 


A LOCAL TECHNIQUE FOR SIMULATED DELAMINATION GROWTH 

In the previous sections, simple quadratic expressions were developed which made it 
possible to calculate the strain energy release rate for various load combinations. In this part of the 
investigation a technique was developed where the forces and displacements at the crack tip (see 
Figure 6) obtained from three linear analyses are superposed. The calculated energy release rates 
for one delamination length are matched with the corresponding results from one nonlinear finite 
element analysis and a correction factor is determined. This correction factor is then used to size the 
results obtained from linear analyses for all other delamination lengths. 

Only one nonlinear finite element analysis was performed for each load case using a full 
model of the damaged specimen as shown in Figure 4. Loads measured at the onset of damage as 
shown in Figure 2 and discussed in the previous paragraph were simulated. Mode I and mode II 
energy release rates G IJJL and G nj ^ were computed for a delamination length equal to the length of 
the first element (a/h =0.181) as shown in Figure 4. Local force and moment resultants N„, Q xy , 
and M were calculated at the location where the end of the frame or stringer flange meets the skin 

xx 

as shown in Figure 11. Resultants plotted in Figure 17 show that the force resultant N xx is zero for 
the three-point bending test as it is free of axial tension. Also as expected, there is a small 
transverse shear, which is non zero. For the tension test, in addition to the membrane resultant, a 
bending moment is present due to the load eccentricity in the flange region and the asymmetric 
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layup of the combined skin and flange laminate with respect to the neutral axis. The shear force 
resultant Q xy is nearly zero, as expected. For the ATB test, calculated membrane and moment 
resultants lie between the computed pure tension and pure bending values [7]. Due to the high 
transverse load during the tests, the shear force resultant is significant for this load condition. It 
was assumed that these local force and moment resultants calculated at the flange tip vary only 
slightly when the delamination is extended. 

Three local sub-models (shown in Figure 18) were then developed to simulate delamination 
growth using a linear analysis. The local sub-model consisted of a small section of the original 
model around the location where the end of the frame or stringer flange meets the skin. To avoid 
any disturbance associated with the load introduction, the length of the model to the left of the 
damage {d x ) was about three times the skin thickness and the length of the model to the right of the 
damage location (<£,) was about three times the skin plus flange thickness (r s +r f ). The mesh used for 
the local sub-model is the same as the mesh of the full model shown in Figure 4. As shown in 
Figure 18(a), boundary conditions for all local sub-models were selected to prevent the translations 
in the plane and rotation of the model. Three unit load cases were simulated as shown in Figures 
18(b) through (d) and the delamination was extended as explained in the paragraph above. External 
loads were chosen such that a unit force resultant Q x or unit moment resultant exists at the 
reference station at the flange tip. For the unit transverse shear load case, a counter reacting 
moment, M c , needs to be applied at the end of the model to assure a pure shear force resultant Q xy 
at the flange tip. To facilitate the simulation of the external moment (Figure 18(c) and (d)) 
three-noded quadratic beam elements with rotational degrees of freedom were used for the 
simulation of the load introduction zone, s, which had die same length as the adjacent plane strain 
elements (Figure 18(a)). A rectangular beam cross section was selected to model the square cross 
section of the skin. The beams were connected to the two-dimensional plane strain model of the 
local section using multi-point constraints to enforce appropriate translations and rotations. This 
procedure was explained for the combined axial tension/bending load case and shown earlier in 
Figure 5. For the beam model, smeared orthotropic material properties were calculated for the skin 
laminate and used as material input data. 

For each unit load case (index N,M,Q), the delaminations were extended and a linear finite 
element analysis was performed for each length a. For each simulation, forces X' Ni (a), X' m (a), 
Xq ( a), and F Ni (a), F Mj (a), F^a), at the delamination tip at node i and the relative displacements 
A^Nnite)’ Au Mm( a )- and Av ' Nm ( a )> Av' Mm (fl), Av’^fat), at the corresponding node m 

behind the delamination tip were retrieved from the finite element results (see Figure 6). Forces at 
node j and relative displacements at node £ were also obtained. In a second step, forces and relative 
displacements for each of unit load cases were scaled by multiplying with the corresponding force 
and moment resultant N^, Q xy and M a obtained from the nonlinear analysis of the full model. The 
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scaled forces and displacements were then superposed yielding 


r,( a )=JV„r N 1 ( a )+M <x -r MI ( a )+e <z r Q 1 (a) 

^ j ( a ) = N„ ' Nj ( a ) "*■ A^xx ’ y Mj ( a ) Qa . ' ^Qj ( a ) 

Av m ( a ) = A^xx ‘ Nm ( a ) + A^xx ' Mm ( a ) Qxz ’ Qm ( a ) (14) 

Av e (a) = N xx • Av (a) + • Av Mi (a) + Q , ^ ■ Av Q/ (a) 


Forces X’j(a) and X' t (a) as well as relative displacements An’ m (a) and A u \(a), were obtained 
accordingly. All forces (X' t (a), X' i (a), and T\{a), Fj(a)), and relative displacements (Au' m (a), 
An \(a), and Av' m (a), Av'/a)) obtained, served as input for the virtual crack closure technique 


Q(a) = - 


2 A a 


y i( a ) ( v m (a)- v ' m -(a)) + y j( fl ) ( v / («)-v>(a)) 


(15) 


AV 'm ( a ) 


AV '| ( a ) 


^ll( a ) — 


2A a 


X\(a)-(u m (a)-u\{a)) + X)(a) (u' t ( a)-u'.(a )) 


Au 'm ( fl ) 


( fl ) 


(16) 


The correction factors c r and c n for mode I and mode n, respectively, were introduced in order to 
size the results for G 1 and G n obtained from the superposition procedure (equations (15) and (16)) 
along the delamination length. One set of correction factors Cj and c n was determined for the entire 
study by matching the G, and G n results obtained for the initial crack (a/h =0.181) with G I>rL and 
Gjjj^l computed from the initial nonlinear analysis. This is- accomplished by calculating G, 
(, a/h =0.181) and G n (a/h =0.181) first with the correction factors set to c r =c n =l and then solving 
for the correction factors 


G tNL (a/ft = 0.181) _O m {a^h z 0181) 

C ’ Gyia/h =0-181) D G n (a//i = 0.181) 

The correction factors obtained for the tension, three-point bending and combined axial 
tension/bending load case are given in Table 2. For the pure tension and the axial tension/bending 
load cases the correction factors are relatively large when compared to the factors calculated for the 
pure bending load case. This is most likely related to the distinct nonlinear load/deflection behavior 
of the specimens subjected to these loadings. Hence, large correction factors are required to match 
the results obtained from the three linear unit load cases with those obtained directly from nonlinear 
FE analysis using VCCT. Consequently, for a nearly linear load/deflection behavior - as observed 
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during the bending test - a much smaller correction factor is required. The load/deformation 
behavior of the specimens for all three load cases is discussed in detail in references 6 and 7. 

For the tension, three-point bending and combined axial tension and bending load case, 
mixed mode energy release rates were calculated using the superposition technique described above 
and given in equations (14) through (17). The results were included in the plots of Figures 14 
through 16. For the initial matrix crack length (a/h =0.181) the results are identical, as this point 
was chosen to match the results and calculate the corrections factors (see equation (17)). The 
correction factors obtained were kept constant during the simulation of delamination growth. The 
obtained mixed mode energy release rates show that G„ increases monotonically for all load cases 
while Gj begins to level off at the longest delamination lengths. For the bending load case the 
results were in excellent agreement with energy release rates calculated directly from nonlinear 
finite element results using VCCT along the entire delamination length. This may be attributed to 
the fact that the load/deflection behavior of the specimen under this load is nearly linear and 
therefore can closely be approximated by the linear analyses of the local sub-models. Along the 
entire delamination length investigated, results were in good agreement for the other load cases as 
well. As the delamination length becomes longer however, the results obtained from the 
superposition technique begin to deviate slightly from the values calculated directly from nonlinear 
finite element analyses. For long delamination lengths it might therefore be advantageous to 
calculate several reference solutions for different delamination lengths from the full model using 
nonlinear analysis and update the corrections factors. 

As mentioned in the previous paragraph, a total of twelve nonlinear analyses were 
necessary when using the conventional approach to obtain the results for one load case as shown in 
Figures 14 through 16. The superposition technique described above required only one nonlinear 
analysis of the full model for each load case and 36 linear analyses of the local sub-model. Even 
for one load case this means a considerable reduction in CPU time. Although additional modeling 
effort is required to create the local sub-model, the results indicate that the proposed technique is 
very efficient for large parametric studies which may occur during preliminary design where 
multiple load combinations must be considered. 


CONCLUDING REMARKS 

Three simple procedures were developed to determine strain energy release rates, G, in 
composite skin/stringer specimens for various combinations of in-plane and out-of-plane loading 
conditions. These procedures may be used for parametric design studies in such a way that only a 
few finite element computations will be necessary for a study of many load combinations. Since 
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energy is a quadratic function of the applied loads, it was not possible to simply superpose and add 
the energy release rates from separate load cases. A simple quadratic expression was previously 
developed to calculate the strain energy release rate for any combination of loads. To validate the 
procedures, results obtained from the quadratic expressions were compared to mode I and 
mode II strain energy release rate contributions, which were calculated from nonlinear two- 
dimensional plane-strain finite element analyses using the virtual crack closure technique. 

For the first technique, the boundary conditions for the tension, bending and combined 
tension/bending load case were chosen in such a manner that the specimen deformation was 
assumed to be a linear function of the applied loads. Therefore a linear finite element solution was 
used to compute the strain energy release rate for various multi-axial load combinations. The 
technique involved solving three unknown parameters needed to determine the energy release rates 
from a simple tension, a simple bending, and one combined tension/bending load case. Excellent 
results were obtained when the external loads were used. This superposition technique, however, 
was only applicable if the structure exhibits a linear load/deflection behavior. 

Consequently, a second modified technique was derived which was applicable also in the 
case of nonlinear load/deformation behavior. The expression derived for the linear case was 
modified such that terms of the external forces were replaced by internal force and moment 
resultants at the flange tip. The energy release rates calculated using this technique seemed 
sufficiently accurate for preliminary design studies. However, force and moment resultants at the 
flange tip need to be calculated and additional effort is required to obtain six unknown parameters 
from a set of six simultaneous linear equations to determine the energy release rates. This 
procedure, therefore, was not time efficient, and hence, less appealing. 

Finally, a third procedure was developed to calculate mixed mode energy release as a 
function of delamination lengths. This procedure required only one nonlinear finite element 
analysis of the specimen with a single delamination length to obtain the force and moment 
resultants at the flange tip and a reference solution for the energy release rates. It was assumed that 
the local force and moment resultants calculated at the flange tip vary only slightly when the 
delamination is extended. Therefore it is sufficient to calculate these resultants for one delamination 
length. The delamination was subsequently extended in three separate linear models of the local 
area in the vicinity of the delamination subjected to unit loads. Forces and displacements computed 
at the delamination tip for the unit load cases were superposed and used in the virtual crack closure 
technique to obtain the distribution of G with delamination length. Results were in good agreement 
with energy release rates calculated directly from nonlinear finite element results using VCCT. 
Although additional modeling effort is required to create the local sub-model, this superposition 
technique is very efficient for large parametric studies which may occur during preliminary design 
where multiple load combinations must be considered. 
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TABLES 



TABLE 1. MATERIAL PROPERTIES. 


E = 1.72 GPa v = 0.30 (assumed isotropic) 


TABLE 2. CORRECTION FACTORS FOR SCALED ENERGY RELEASE RATES. 


Tension Load Case 

Bending Load Case 

Axial Tension/Bending 
Load Case 

ci =1.2657 

ci =1.0036 

ci =1.2791 

c n =1.2484 

cn =1.0646 

c n =1.1720 
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Figure 1. Specimen Configuration. 
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Corner 4 


Comer 2 



Comer 3 Corner 1 

(a) Specimen with Crack Locations. 



(b) Comers 1 and 4 


Delamination B1 Delamination B2 



(c) Comers 2 and 3 


Figure 3. Typical Damage Patterns [6,7] 
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specimen modeled with top grip, axial load cell and pin modeled 

2D plane strain elements with beam elements (E=210 GPa, v =0.3) 




2 nodes with identical coordinates 
beam node A 
2D quad node 15 

multi-point constraints: 

Ua = Uis, v A = v 15 
4>a = ( u a - u, ) /t s 
u, = u, + y, ( Uzg - u, ) / U 
at x =1 01 .6 mm 


Figure 5. Loads and Boundary Conditions for the 
Combined Axial Tension and Bending Test. 
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P= 5.5 kN 
11.0 kN 
16.5 kN 

x,u,P 22.0 kN 

(a) Tension Load Case 


Q=1 12.5 N 337.5 N 



(b) Bending Load Case 




x.u.P 


P= 5.5 kN 

1 1 .0 kN 
16.5 kN 

22.0 kN 


(c) Combined Load Condition 


Figure 7. Loads and Boundary Conditions For Tension and Three- 
Point Bending and Combined Loading Case. 
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Figure 8. Comparison of Computed Strain Energy Release Rates 
with Superposed Values for Tension Load Case. 
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Figure 9. Comparison of Computed Strain Energy Release Rate Components with 
Superposed Values for Combined Tension and Bending Load Case. 
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Applied Transverse Load Q, N 


Figure 10. Comparison of Computed Total Strain Energy Release Rates with 
Superposed Values for Combined Tension and Bending Load Cases. 



Figure 1 1. Calculation of Force and Moment Resultants 
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0 5 10 15 20 25 30 35 


Applied Transverse Displacement v, mm 

Figure 12. Comparison of Computed Strain Energy Release Rate Components 
with Scaled Values for Combined Tension and Bending Load Case. 
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Applied Transverse Displacement v, mm 

Figure 13. Comparison of Computed Total Strain Energy Release Rates with 
Scaled Values for Combined Tension and Bending Load Cases. 
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Delamination Length a, mm 


Figure 14. Computed Strain Energy Release Rates for Delamination Growth 
in a 90745° Flange Ply Interface for Tension Load Case. 



Figure 15. Computed Strain Energy Release Rates for Delamination Growth 
in a 90745° Flange Ply Interface for Three-Point Bending Load Case. 
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Delamination Length a, mm 

Figure 16. Computed Strain Energy Release Rates for Delamination Growth in 
a 90°/45° Flange Ply Interface for Combined Tension and Bending Load Case. 



N mm/mm 



Figure 17. Computed Force and Moment Resultants at Flange Tip. 
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obtain forces 
YV. 1 Y , „ 1 XV..XV, 
and displacements 

V'Nm, V Nm . , V'w, V' N /• 
UNm, U' Nm ., U' W , U'n/* 
for VCCT 

(b) Unit axial tension load case 


M,<|> 

r 

y.v.Q 


► x.u.N 


(a) Local finite element model 




obtain forces 
Y Ml , Y'mj , X'mi , X*Mj 
and displacements 
Ymitii Y Mm* , V'm<, V'mT 
UMm, U'Mm- . U'm/, U'mT 

for VCCT 


obtain forces 
Yqi, Y Qi ,X' Q ,X' Qj 
and displacements 

Yom, Yom*, V'q {, V‘q f 
Ll'orm U'om* , Uq/, U'q/* 

for VCCT 


(d) Unit tranverse shear load case 

Figure 18. Local finite element model for linear analyses and unit loads. 
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